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We develop a numerical procedure to efficiently model the nonequilibrium steady state of one¬ 
dimensional arrays of open quantum systems, based on a matrix-product operator ansatz for the 
density matrix. The procedure searches for the null eigenvalue of the Liouvillian superoperator 
by sweeping along the system while carrying out a partial diagonalization of the single-site sta¬ 
tionary problem. It bears full analogy to the density-matrix renormalization group approach to 
the ground state of isolated systems, and its numerical complexity scales as a power law with the 
bond dimension. The method brings considerable advantage when compared to the integration of 
the time-dependent problem via Trotter decomposition, as it can address arbitrarily long-ranged 
couplings. Additionally, it ensures numerical stability in the case of weakly dissipative systems 
thanks to a slow tuning of the dissipation rates along the sweeps. We have tested the method on a 
driven-dissipative spin chain, under various assumptions for the Hamiltonian, drive, and dissipation 
parameters, and compared the results to those obtained both by Trotter dynamics and Monte-Carlo 
wave function. Accurate and numerically stable convergence was always achieved when applying 
the method to systems with a gapped Liouvillian and a non-degenerate steady-state. 
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I. INTRODUCTION 

The study of the nonequilibrium dynamics of open 
many-body quantum systems has gained significant mo¬ 
mentum in recent years, thanks to the experimental 
progress achieved in several areas, including ultracold 
atoms in optical lattices m, trapped ions [SHiQ] , arrays 
of optical micro-resonators and superconducting 

circuits pMUTB] , 

A feature common to all these systems is the coupling 
to an external environment in the form of coherent or in¬ 
coherent input and output channels. The time evolution 
of the system is then governed by an interplay of the 
Hamiltonian and the driven-dissipative dynamics. For 
stationary external conditions, this dynamics typically 
leads to a nonequilibrium steady state (NESS), for which 
a multitude of novel phenomena are expected, including 
nonequilibrium quantum phase transitions UMB and 
the possibility of engineering quantum states through 
tailored dissipation [22], in view of advanced quantum 
information strategies [23] . 

The theoretical description and modeling of open 
quantum systems out of equilibrium represents a ma¬ 
jor challenge. Indeed, similarly to the ground state of 
isolated many-body quantum systems, the NESS can be 
characterized by quantum correlations which - partic¬ 
ularly when approaching criticality - require for their 
exact determination a computational effort that scales 
exponentially with the system size [24l [25]. As added 
difficulties however, the NESS is generally not a pure 
quantum state, nor can it be directly determined from 
the Gibbs principle, as in the case of thermal equilib¬ 
rium. 

Generally, an open quantum system is described by a 
density matrix p, whose dynamics obeys the Von Neu¬ 
mann equation p = Cp dictated by the Liouvillian su¬ 


peroperator C (we set h = 1 here and in what fol¬ 
lows) [26| [27] . Two strategies are then available for the 
determination of the NESS. Eirst, one can directly in¬ 
tegrate the time evolution until stationarity is reached. 
Second, a solution of the equation Cp = 0 can be directly 
computed, under the additional condition that Tr(p) = 1. 
Apart from special cases in which analytical solutions can 
be found [28l [29] , both strategies can be handled numer¬ 
ically only for very small systems [3QH3^ - if an exact 
solution is sought. Larger systems typically require some 
level of approximation and, still in recent times, many 
studies have restricted to mean-field approximations [34]- 
[38] . thus neglecting quantum correlations. Only very re¬ 
cently a variational principle for the NESS of open quan¬ 
tum systems has been demonstrated [39] and applied to 
1-D systems 001, while a spatial decimation method spe¬ 
cific to the stationary Von Neumann problem has been 
proposed [4T] . 

In this scenario, one-dimensional systems represent 
a special case in which a very accurate description of 
the many-body quantum state is made possible thanks 
to the advent of the Density Matrix Renormalization 
Group [42H^ (DMRG) and of the equivalent variational 
approach based on the Matrix Product State (MPS) 
ansatz [45] [46]. In typical situations, the MPS-DMRG 
approach allows a surprisingly good account of quantum 
correlations at finite spatial range, with a computational 
overhead that scales polynomially with the dimension of 
the Hilbert space. The MPS approach has been suc¬ 
cessfully extended to the modelling of the unitary time 
evolution of a closed quantum system [45l [46] . Eor open 
quantum systems, an analogous Matrix Product Opera¬ 
tor (MPO) ansatz for the density matrix has been pro¬ 
posed and applied to model both thermal equilibrium [47] 
and temporal dynamics [471[48]. In particular, the long 
time dynamics has been employed to obtain the nonequi- 
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librium steady state (NESS) - in presence of driving fields 
and dissipation - in different settings [49l [54H58] . There 
are however several settings in which the MPO dynami¬ 
cal approach to the NESS suffers from limitations. This 
is the case in presence of slow dissipation rates (com¬ 
pared to the energy scale set by the Hamiltonian), or 
when dissipation acts on a small part of the system only, 
as in transport configurations [sniEaisn]- Some systems 
may even display an algebraic, rather then exponential 
dynamics to the NESS j49f[53] . Einally the Trotter de¬ 
composition, typically used in these dynamical schemes, 
suffers from a severe limitation: it is restricted to nearest 
neighbor couplings. Only recently, numerical approaches 
have been suggested [6QH63] . that overcome this limita¬ 
tion, but only for modeling the unitary dynamics of iso¬ 
lated systems. 

Here, we develop an efficient implementation of the 
variational principle to directly determine the NESS of 
Markovian open quantum systems. The method does not 
rely on the integration of the long-time dynamics, thus 
lifting all the limitations described above. The varia¬ 
tional principle for determining the NESS has recently 
been proposed isiiioi and implemented within an MPS- 
DMRG scheme [40] . The approach that we propose relies 
directly on the search for the zero-eigenvalue of the su¬ 
peroperator C for the determination of the NESS. This 
approach has shown full numerical stability when applied 
to gapped Liouvillians with a non-degenerate NESS. As a 
test of the method, we simulate a driven-dissipative Ising 
chain, and compare the results to those obtained by simu¬ 
lating the MPO dynamics [54l|64| and with Monte-Carlo 
Wave Eunction (MCWE) [651468] . We then simulate the 
same system, in presence of longer-range couplings or 
slow dissipation rates, thus showing its wide range of 
applicability in the description of driven dissipative sys¬ 
tems. We finally discuss the computational complexity 
of the approach and compare it to other existing meth¬ 
ods [40]. 

II. THE METHOD 

We consider a one-dimensional chain of N coupled 
quantum systems, each characterized by d possible states, 
in the presence of external driving fields and Markovian 
coupling to the external environment. The dynamics is 
governed by the Lindblad-Von Neumann master equa¬ 
tion [2i [27] 

^ = Cp = -i[H, p]-IY1 , 

i 

( 1 ) 

where Ki are the operators corresponding to the transi¬ 
tions induced by the environment. The NESS solution 
obeys the equation >Cpness = 0- 

Eor the purpose of numerical implementation, it is con¬ 
venient to map the MPO representation onto an equiv¬ 
alent MPS form. We do this by the vectorization pro¬ 


cedure, where the density matrix p is reshaped into a 
column vector, here denoted by |p)), by concatenating 
all its columns. To express the Liouvillian superop¬ 
erator in this representation, we rely on the property 
\XpY)) = (8) X|p)), where X and Y are matrices. 

Then, C takes the form of the matrix defined by [69] : 

c = -i(I ®1) 

+ i ®ki-i®klki-kfk* ®i). (2) 

The determination of the NESS can then be reformu¬ 
lated as the variational minimization of the euclidean 
norm functional 

ll>C|/3))||>0. (3) 

The MPO representation of the density matrix reads 

I(Ji... (7/... (Tat) are the states of the system, \aj) is the 
state of the j-th site of the chain and the sets of matri¬ 
ces {A} parametrizes the MPO state [45j|46]. Through 
vectorization, we may express the density matrix as an 
MPS 


where |5])) = ||cr)(cr'|)), and the indices in the matrix 
elements of p have been encoded as = (a- — l)d + cr^. 
Once expressed using a MPS representation, the problem 
is determined by @ can be solved using the MPS-DMRG 
strategy, for which we will refer to the treatment - and 
the related notation - extensively presented in Refl46l In 
particular, in order to derive the equation for the on-site 
problem, it is useful to express the density matrix in a 
mixed canonical form [46] 

1-5))=E p""' (n 1^)) ’ (5) 

s \i=l / \i=Z+l / 

where the matrices , associated to the l-th site in the 
MPS ansatz, have been singled out from the MPS expres¬ 
sion, and the sets of matrices {A} and {H}, with maximal 
bond dimension D, are left and right normalized, respec¬ 
tively [45[ [46]. The MPS is depicted in Eig. [^a) in the 
usual diagrammatic representation gSlISD]- The symbol 
g then denotes a rank-three tensor, and is associated to 
a local representation of the density matrix at site 1. 

The Liouvillian operator can be represented in an 
MPO form as 

N 

E (6) 

i=i 

as depicted in Eig.j^b). Here, Dw is the bond dimension 
of the MPO representation of £, i.e. the dimension of the 
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FIG. 1. Diagrammatic representations of matrix products, 
(a) Diagram of the vectorized density matrix as an MPS in 
the mixed canonical form 0. (b) Diagram of the Liouvillian 
operator in the MPO representation (c) Diagram repre¬ 
senting the on-site Liouvillian operator Hi. In all diagrams, 
triangles pointing right represent the left-normalized matrices 
A, while triangles pointing left denote the right-normalized 
matrices B, both entering the mixed canonical form of the 
MPS, Eq. if- The local representation of the density matrix 
Q is depicted as a circle at site 1. Thin lines represent physical 
indices while thick lines denote bond indices. The MPO ma¬ 
trices W in Eq. ^ are represented as squares, (d) Diagram¬ 
matic scheme illustrating the computational complexity asso¬ 
ciated with index contractions at one site, in the case where 
the simple Liouvillian C is used in the variational approach, 
(e) Same as (d), in the case where the squared Liouvillian C 
is instead used. In this second case, contraction of the MPO 
bond indices bears an additional 0{Dw) computational cost. 


matrices W in Dw is defined by the complexity of 
the system Hamiltonian and dissipative processes and is 
fixed for a given model [46] . 

A most natural choice for the variational determi¬ 
nation of the NESS, as adopted in Ref. [40l would be 
to express ^ as {{p\C^C\p)). In this way, the prob¬ 
lem bears a full analogy to the MPS-DMRG approach 
to isolated systems, with the hermitian, semi-positive- 
defined operator C playing the role of the Hamilto¬ 
nian. However, this choice requires handling the prod¬ 


uct jC^jC at some level within the algorithm. Let us as¬ 
sume a given MPO representation of £, with bond di¬ 
mension Dw‘ As depicted in Fig. Bd), when computing 
the quantity {{p\C\p)) the numerical complexity associ¬ 
ated to the index contractions on each site scales with 
The corresponding complexity, in the case of 
the quantity {{p\C^C\p)), is sketched in Fig. [^e) and 
would naively scale as 0{D^). More specifically, the 
computational complexity according to Ref. [46] (see Eq. 
(197)) is 0{2dPD^Dw + for the C algorithm 

and would be 0{2(PD^D‘^ -^d‘^D‘^D^) for the C^C case 
if one opted for directly using a MPO representation of 
the squared Liouvillian of bond dimension D‘^. However, 
rather then constructing the C MPO, one may improve 
the second approach by storing only C and carrying out 
the matrix multiplication by “on the fly” at each op¬ 
timization step. This would reduce the complexity of the 
jC^jC approach to 0{2{(PD^D‘^^d^D‘^D^), which is still 
however one Dw-^SiCtoi slower than the C approach. For 
models such as bilinear-biquadratic hamiltonians, or even 
XYZ models with slightly involved dissipative processes, 
the MPO representation of C can reach bond dimension 
easily exceeding Dw ~ 10. The present strategy may 
thus easily lead to a computational gain of more than 
one order of magnitude. Furthermore, in most systems 
of interest, £ is a very sparse matrix, and this computa¬ 
tional advantage is partly spoiled when instead adopting 
the generally less sparse squared Liouvillian. Finally, the 
bond dimension of the Liouvillian MPO has a relevant 
computational impact also on the iterative solution of 
the on-site eigenvalue problem at each site of the chain, 
in cases where the matrix is not fully stored and the lin¬ 
ear operator is instead applied to vectors in a functional 
fashion. In these cases, for an MPO with bond dimension 
Dw^ the complexity associated to the matrix-to-vector 
multiplication is 0{D^Dwd^ ^D^Dl^d^) (see 

equation (201) of [46]), again highlighting the importance 
of using an MPO with minimal bond dimension. These 
considerations led us to explore the possibility of finding 
the NESS by directly searching for the null eigenvalue of 
C. 

In the MPS-DMRG algorithm, all matrices A and B in 
(§ are kept constant, and the inequality (§ can then be 
cast into an on-site linear problem for the optimization 
of g. To this purpose, we introduce the on-site Liou¬ 
villian operator Ci for site /, which is a rank-six tensor 
obtained from the quantities C and \p)) by contracting all 
indices associated to the other lattice sites, as depicted in 
Fig. I^c). The minimization of the norm functional § 
is then achieved by solving the local problem Ciq = ^ 
successively on each site of the chain, sweeping along 
the chain in both directions until convergence to the null 
eigenvalue is reached. To this purpose, we solve the lo¬ 
cal problem by computing the complex eigenvalue of Ci 
closest to a small target scalar value. This scalar must 
be chosen much smaller than all energy scales character¬ 
izing the problem, in order to achieve convergence to the 
null eigenvalue of C. Convergence is achieved after a suffi- 
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cient number of sweeps, and by choosing bond dimensions 
large enough to accurately model the quantum correla¬ 
tions arising in the NESS. For the eigenvalue problem, 
we adopted here the Shift-and-Invert Arnold! method, 
which is most efficient for small magnitude eigenvalues. 
The method has yielded in our tests the most stable and 
efficient realization of the algorithm. Due to the matrix 
inversion however, the Shift-and-Invert method requires 
the full storage of the local Liouvillian, i.e. a memory cost 
0{D‘^(P X D‘^(P). In cases where this memory cost can¬ 
not be afforded, it is still possible to adopt direct iterative 
schemes. The ARPACK library in particular m makes 
non-inverting versions of the Arnold! method available. 
Our experience is that, while solving the storage problem, 
these methods are generally considerably slower and less 
stable - though only marginally - than the shift-and- 
invert method. 

In general, a matrix diagonalization targeting small 
complex eigenvalues is usually characterized by slow con¬ 
vergence. To overcome this limitation in the present case, 
and increase efficiency, we start the computation using a 
small bond dimension, and allow it to increase gradually 
along the sweeps, by each time padding the larger den¬ 
sity matrix with zeros. Lastly, we have found that the 
algorithm could become unstable when directly targeting 
very small dissipation rates (compared to the Hamilto¬ 
nian energy scale). To ensure the stability of our imple¬ 
mentation in such cases, we start the computation using 
larger dissipation rates, and let them decrease exponen¬ 
tially towards the desired values along the sweeps. In 
practice, in our tests we started the computation with 
values of D between 5 and 10 and run several tens of 
sweeps, while gradually decreasing the dissipation rates 
if needed. We observed that this first phase can be sped 
up significantly by restricting the number of iterations of 
the Shift-and-Invert algorithm to less than 10. After this 
first phase has converged, we refine the result by allow¬ 
ing the bond dimension to increase gradually, while at 
the same time increasing the number of Shift-and-Invert 
iterations in each step to a few hundreds. This second 
phase typically requires less than 10 sweeps to achieve 
full convergence. 

Since the introduction of MPS modelling of mixed 
states, the issue of preserving the positivity of the density 
matrix has been discussed m and shown to be NP-hard 
to verify E). It should be noted that only very recently 
a local purification scheme for the Trotter evolution has 
been proposed in which guarantees positivity of the 
density matrix. However, we stress that for the cases 
we have considered we never encountered convergence to 
an MPS that presented unphysical results and we never 
had to reinforce the density matrix properties which, in¬ 
stead systematically result from the convergence of the 
algorithm. 

Note also, that in the MPS approach the state is nor¬ 
malized according to the euclidean norm, i.e. {{p\p)) = 1 . 
Thus, in general, the condition on the trace Tr(p) = 
((IIp)) = 1 is not automatically fulfilled, and the expecta¬ 



FIG. 2 . Comparison between the spatial correlations 
(Xm^m+i) for I = 1,...,4, computed by Trotter dynamics 
and through the direct variational MPS determination of the 
NESS. Parameters were set as E = 0 and 7 = J, and the 
array length was A = 15. For the Trotter dynamics, the time 
step was set as dt = O.IJ and time integration was carried 
out until T = 10 / 7 . For both the Trotter dynamics and the 
variational NESS methods the bond dimension was at most 
D = 20. These results also displayed perfect agreement with 
MCWF calculations (not shown). 


tion value of an arbitrary observable O must be evaluated 
as ( 6 ) = Tr(pO)/Tr(p) = ((I|I ^ 6 |p))/((I|p)). 


III. RESULTS 


As a test of the method, we simulate a driven- 
dissipative quantum Ising chain [SH [64], described by 
the Hamiltonian 

^ = E[ hZi + JXiXi+i + VXiXi+2 

i 

with h being a local effective magnetic field, and J and 
V respectively the coupling between nearest neighbours 
and next nearest neighbours. The dissipative part is pro¬ 
vided by transition operators K = — iY){2 at each 

site, with X, Y and Z being Pauli matrices and 7 the dis¬ 
sipation strength. 

Our study focuses on three paradigmatic cases and we 
initially assume a small system size (15 sites), to allow 
for a direct comparison with MCWF simulations. The 
MCWF unravels the master equation for density ma¬ 
trix into stochastic pure state trajectories in the Hilbert 
space. Dissipation is accounted for by non-hermitian 
terms in the Hamiltonian, while the corresponding fluctu¬ 
ations are enforced by random “quantum jumps” gener¬ 
ated with a probability proportional to the square root of 
each dissipation rate. The method is described in detail 
in [65j|66], and we specifically adopted the QuTiP tool¬ 
box El for all MCWF calculations. We also compare to 
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FIG. 3. Comparison between the spatial correlations 
(XrnXm+i) for / = 1,..., 4, computcd by MCWF and through 
the direct variational MPS determination of the NESS. Pa¬ 
rameters were set as V = 0.5J and 7 = J, and the array 
length was N" = 15. The MCWF simulation was performed 
with 1000 trajectories and time integration was carried out 
until T = 10 / 7 . For the variational NESS calculations, the 
bond dimension was at most D = 30. 


the standard Trotter MPS evolution [54] for benchmark¬ 
ing the method. 

In the first case, nearest neighbour couplings are con¬ 
sidered, as in Ref. [ 54 J In Fig. Q, the results for the cor¬ 
relations {XrnXm-\-i) I = 1,..., 4, obtained both using 
Trotter dynamics and the variational method are shown. 
They coincide perfectly with each other and with the 
data obtained in Ref. [54j In particular, the system dis¬ 
plays ferromagnetic order for negative external field and 
anti-ferromagnetic order for positive external field. The 
small discrepancy observed between the data obtained 
with the two methods is simply due to the Trotter error. 
For this case, 7 = J and the driven-dissipative time evo¬ 
lution is well handled by the Trotter dynamics, which is 
therefore the method of choice, as the time scale to reach 
the NESS is short and the resulting simulation turns out 
to be much faster than the variational method. This con¬ 
sideration holds in general, in cases with next-neighbour 
couplings and sufficiently fast dissipation rates. 

The second case we study, is that of a system with 
longer range couplings. In this case the usual Trotter 
dynamics cannot be employed and thus the variational 
NESS becomes the natural method of choice. In Fig. 
we compare results obtained with the MCWF and vari¬ 
ational methods. Once again, we obtain a very good 
agreement between the two methods, even for small bond 
dimension. The next nearest neighbour coupling ampli¬ 
fies the ferromagnetic correlations, while having a size¬ 
able effect on the anti-ferromagnetic side. By comparing 
Fig-i and Fig. [^we see that the next-nearest-neighbour 
correlation (X^X^+ 2 ) changes sign and we only observe 
anti-correlation at longer distance ((X^X^+ 3 )). We ar¬ 
gue that, when adding genuinely long ranged couplings. 



FIG. 4. Comparison between the spatial correlations 
(XrnXm+i) for / = 1,..., 4, computcd by MCWF and through 
the direct variational MPS determination of the NESS, in the 
case of weak dissipation. Parameters were set as 7 = O.IJ, 
and the array length was N = 15. The MCWE simulation was 
performed with 16 trajectories and time averages were taken 
from t = 10/7 and until t = IOO/ 7 . Eor the variational NESS 
the bond dimension was at most D = 50. The inset shows 
the correlations {XrnXm+i) as a function of I with h = J 
for a system of N" = 50 sites and bond dimension D = 60. 
This system size lies beyond the computational reach of the 
MCWE method. 

the anti-ferromagnetic order in the positive external field 
sector might be completely suppressed. 

As the third case, we simulate the same model in pres¬ 
ence of a small dissipation rate. In this case, dynamical 
methods will become less effective and converge slowly. 
We have observed that the variational method in this 
case could become unstable. This issue was however 
completely removed by adopting a gradual decrease of 
the dissipation rate along the sweeps, as discussed pre¬ 
viously. In this case, the small dissipation rate results 
in increased correlations, both in the ferromagnetic and 
anti-ferromagnetic case, as show in Fig.[^ It is also inter¬ 
esting that nontrivial correlations emerge for very small 
external field showing that there are still novel regimes 
to be explored for these driven dissipative systems. The 
inset in Fig. shows the correlations {XrnXm-\-i) as a 
function of /, computed for a longer system with = 50 
sites. The combination of a quasi-local hamiltonian with 
an on-site dissipation mechanism seems to generally lead 
to an exponential decay of the correlations. This setting 
typically holds for driven dissipative optical systems such 
as coupled optical cavities. This result suggests that the 
present method may efficiently model the NESS of long 
one-dimensional systems, already at moderate bond di¬ 
mension. 


IV. CONCLUSION 

In conclusion, we have presented an efficient imple¬ 
mentation of the variational principle for the NESS of 
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one-dimensional driven-dissipative quantum systems us¬ 
ing an MPO ansatz for the density matrix. The compu¬ 
tational overhead of the method scales as a power law 
both in the dimension of the Hilbert space and in the 
bond dimension of the MPO. Vectorization allows to map 
the problem onto an effective linear eigenvalue problem, 
that can be then solved using a MPS-DMRG approach. 
We have applied the method to a model spin chain as a 
test, under various assumptions for the parameters. As 
compared to direct integration of the system dynamics, 
the present approach brings considerable advantage in 
cases where the dissipation rates are slow compared to 
the Hamiltonian energy scale. In particular, through a 
slow tuning of both the MPS bond dimension and the dis¬ 
sipation rates towards the target values, numerical sta¬ 
bility and convergence to the physical NESS is achieved 
in all cases that we have studied. Also, the method gives 
access to systems with long-range couplings, for which 
the standard Trotter dynamics cannot be employed. In 
such cases, new algorithmic approaches to direct time in¬ 
tegration have very recently emerged [sniini]. The direct 
comparison between the present approach and these new 
developments is left as a venue for future investigations. 

Modeling nonlinear driven-dissipative quantum sys¬ 
tems generally represents a major challenge, as these sys¬ 


tems combine the inherent difficulty in correctly describ¬ 
ing quantum correlations to the nonequilibrium character 
of their approach to stationarity. This difficulty emerges, 
in particular, when dynamical critical phenomena and 
quantum phase transitions occur. Then, quantum cor¬ 
relations typically acquire a long spatial range and may 
even decay algebraically m- Methods relying on the 
MPS ansatz are in these cases an ideal tool, as they 
provide control over the spatial range of quantum cor¬ 
relations through the bond dimension, while preserving 
a power-law computational complexity. In this frame¬ 
work, the method presented in this work holds promise 
as a powerful tool for the study of emergent quantum 
phenomena in nonequilibrium open quantum systems. 
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